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Abstract: High Mountain Asia (HMA) region contains the world's highest peaks and the largest 
concentration of glaciers except for the polar regions, making it sensitive to global climate change. In the 
context of global warming, most glaciers in the HMA show various degrees of negative mass balance, 
while some show positive or near-neutral balance. Many studies have reported that spatial heterogeneity in 
glacier mass balance is strongly related to a combination of climate parameters. However, this spatial 
heterogeneity may vary according to the dynamic patterns of climate change at regional or continental 
scale. The reasons for this may be related to non-climatic factors. To understand the mechanisms by which 
spatial heterogeneity forms, it is necessary to establish the relationships between glacier mass balance and 
environmental factors related to topography and morphology. In this study, climate, topography, 
morphology, and other environmental factors are investigated. Geodetector and linear regression analysis 
were used to explore the driving factors of spatial variability of glacier mass balance in the HMA by using 
elevation change data during 2000—2016. The results show that the coverage of supraglacial debris is an 
essential factor affecting the spatial heterogeneity of glacier mass balance, followed by climatic factors and 
topographic factors, especially the median elevation and slope in the HMA. There ate some differences 
among mountain regions and the explanatory power of climatic factors on the spatial differentiation of 
glacier mass balance in each mountain region is weak, indicating that climatic background of each 
mountain region is similar. Therefore, under similar climatic backgrounds, the median elevation and slope 
are most correlated with glacier mass balance. The interaction of various factors is enhanced, but no 
unified interaction factor plays a primary role. Topographic and morphological factors also control the 
spatial heterogeneity of glacier mass balance by influencing its sensitivity to climate change. In conclusion, 
geodetector method provides an objective framework for revealing the factors controlling glacier mass 
balance. 
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1 Introduction 


Temporal and spatial variation characteristics of glacier mass balance are closely related to the 
climate and environment. Glacier mass balance affects a series of glacier material properties (such 
as ice formation, temperature and movement), and glacier scale (such as area, length, and volume). 
Therefore, glacier mass balance is the critical link to the climate environment. High Mountain 
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Asia (HMA) is the most glacially developed region in the middle latitudes and the largest glacier 
concentration except for the polar regions (Farinotti, 2017). As a result of interaction with 
westerlies, monsoons, and global warming, glacier mass changes will significantly impact 
downstream water resources and water environment (Yao et al., 2012). However, the change of 
glacier mass balance is spatially heterogeneous (Yao et al., 2012; Brun et al., 2017; Hugonnet et 
al., 2021). Most glaciers show negative mass balance, while a few regions (such as the West 
Kunlun and Karakoram mountains) maintain stable or positive mass balance (Brun et al., 2017). 
The reasons for these temporal and spatial variations remain unclear, as topography and other 
environmental factors complicate the glacier's response to climate change (Brun et al., 2019). 

There are two possible explanations for the differences in glacier mass balance: spatial 
heterogeneity of climate change and different glacier responses to climate change (Sakai and 
Fujita, 2017). Recent studies have shown that glacier mass in the Nyaingéntanglha and Hengduan 
mountains has decreased dramatically, while glacier mass in the Karakoram and West Kunlun 
mountains has increased slightly (Brun et al., 2017). The latter increase, known as the 
"Karakoram Anomaly" (Hewitt, 2005), is thought to be due to cooling temperatures in summer, 
and increased precipitation in winter (Farinotti et al., 2020). However, the effects of climate 
heterogeneity across the HMA have not been fully quantified, leaving the possibility that most of 
anomalies are due to the glaciers' different responses to similar climate changes. A simple index, 
such as the sensitivity of glacier mass balance (at the individual glacier or regional scale) to 
temperature and precipitation (Oerlemans and Fortuin, 1992; Marzeion et al., 2018) can be 
numerically estimated heterogeneity in the response of glacier mass balance to climate change. 
Sakai and Fujita (2017) suggested that, in explaining the recent glacial changes in the HMA, 
spatial heterogeneity in glacier mass balance is more informative than regional difference in 
climate change. Whether or not the local or regional sensitivity of glacier mass balance to 
temperature and precipitation is related to the topographic parameters of glacier is not well 
understood. 

Furthermore, the sensitivity of glacier to climate change depends on the current environment of 
the glacier, i.e., the environmental control over the heat and mass balance of the glacier. Salerno 
et al. (2017) conducted a statistical analysis of the relationship between glacier thinning rates and 
morphological variables in the Himalaya Mountains. They found that reduced downstream 
surface slope was the primary morphological variable controlling the thickness changes of 
glaciers. Scherler et al. (2011) demonstrated that supraglacial debris had caused spatial variations 
of the glacier fluctuations in the Karakoram and Himalaya mountains. Brun et al. (2019) reported 
that the glacier tongue slope, mean elevation, debris cover, and the avalanche contribution area 
together account for the mass balance variability in the HMA. Few studies have investigated the 
controlling factors of glacier mass balance in the HMA, and most are focused on local areas. This 
kind of research for other areas outside the HMA also mainly focuses on local area research. In 
the European Alps Mountains, the mass balance of an individual glacier is also related to glacier 
morphology (Huss, 2012; Fischer et al., 2015). Potentially, the outcomes from previous studies 
investigating the controlling factors of glacier mass balance are sometimes contradictory, just 
because the factors controlling the spatial distribution of glacier mass balance vary from region to 
region (i.e., the outcomes and findings from one region cannot be generalized, and applied to 
other regions). Moreover, previous studies used linear regression, making it difficult to detect any 
nonlinear relationships between environmental variables and glacier mass balance. The 
relationships between glacier mass balance and environmental variables are often not simply 
linear due to the complex dynamic processes involved (Li et al., 2019). Determining such 
relationships, especially the potentially nonlinear ones, is helpful in understanding the influences 
of environmental factors on glacier mass balance. 

Geodetector can detect the linear and nonlinear relationships between glacier mass balance and 
environmental variables (Wang et al., 2010). The current study uses geodetector method to 
quantitatively analyze the explanatory power of climate, topography, and other environmental 


ZHANG Zhen et al.: Spatial variability between glacier mass balance and... 


factors on glacier mass balance. These relationships are considered not only from a whole HMA 
perspective but also by investigating individual mountain region. Overall, the driving mechanism of 
the spatial heterogeneity of glacier mass balance in the HMA is examined from a new perspective. 


2 Data and methods 


2.1 Study area 


HMA refers to the high-altitude region of Asia composed of the Tibetan Plateau and its 
surrounding areas. It spans several Central Asian countries within 23°N-47°N, 65°E-105°E, from 
the Pamir Plateau and Hindu Kush Mountains in the west to the Hengduan Mountains in the east, 
the Tianshan Mountains in the north, and the Himalaya Mountains in the south. HMA is the 
highest mountain area on Earth, and is the largest glacier-gathering area except for the polar 
regions. It is also the largest glacial region of the mid-to-low latitudes (Farinotti, 2017). Known as 
the "third pole" of the world, the area is sensitive to global climate change, and has a vulnerable 
ecosystem with an extensive area of glaciers, permafrost, snow, lakes, and land. The ecosystem 
makes HMA a unique geographical unit; it not only regulates runoff from glacial meltwater but 
also influences climate change at the Asian and even global scales, thereby being an indicator and 
regulator of global climate change (Yao et al., 2015). Therefore, HMA plays very important roles 
in ecological security, human survival, and sustainable economic development in Asia. 


2.2 Glacier mass balance data 


Mass balance data for individual glacier come from regional statistics of ArcGIS software-based 
on geodetic mass balance, which was extracted from glacier elevation change rate during 
2000-2016 as described by Brun et al. (2017). Only glaciers with an area larger than 2 km? were 
included in this study to minimize errors. We used Randolph Glacier Inventory v. 6.0 (RGI 6.0) as 
glacier mask data. 


2.3 Glacier morphological data 


The area (A), minimum (Zmin), and maximum elevation (Zmax), median elevation (Zmea), slope (S), 
aspect (As), and length (Lmax) of glaciers were derived from the attribute data of RGI 6.0. The 
elevation difference between Zmax and Zmin is called the difference of glaciation (DG; i.e., Zmax-Zmin). 
The elevation difference between Zmax and snowline altitude of a glacier (here taken as Zmea) is 
called the positive difference of glaciation (PDG; i.e., Zmax-Zmea). Corresponding to PDG, the 
negative difference of glaciation (NDG) is defined as Zmea-Zmin. Hypsometric index (HI) is 
calculated by Equation 1. 


A = 
HI = óm Ame 56. < HI <1, then: I=, (1) 
Z HI 


med min 
The debris cover dataset (Scherler et al., 2018) provides the global distribution of supraglacial 
debris cover, and is available on GFZ (Helmholtz Centre Potsdam, German Research Centre for 
Geosciences) Data Services (https://dataservices.gfz-potsdam.de/panmetaworks/showshort.php?id= 
escidoc:3473902). The debris coverage rate (DCR) is the proportion of supraglacial debris area to 
the total glacier area. 


2.4 Meteorological data 


We used monthly ERAS data (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land- 
monthly-means?tab=overview) to calculate the annual mean temperature (7), mean annual 
precipitation (P), the changing trend of annual mean temperature (T.), and mean annual 
precipitation (P.) from 2000 to 2016. The changing trend of T and P was derived using linear 
regression. Since the resolution of ERA5 data is 0.1°, most glaciers are smaller than a grid area. The 
meteorological data of glaciers smaller than a grid area are taken from the grid value, and values of 
glaciers larger than a grid area are calculated using the regional statistics method in ArcGIS 
software. 
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2.5 Geodetector method 


Previous studies utilized spatial regression methods (such as correlation analysis), which linearly 
correlated dependent variable Y to independent variable X. The heterogeneity of the spatial 
distributions of dependent variable Y, and independent variable X also reflects the correlation 
between the two variables, which includes both linear and nonlinear portions. Both correlation 
relationships can be included and considered within geodetector. When linear regression is 
significant, geodetector must be significant. When linear regression is not significant, geodetector 
may still be significant as long as there is a relationship between the two variables that geodetector 
can characterize. The detailed information on geodetector can be referenced from the research of 
Wang et al. (2010). Following is a brief review of the application of this approach to the current 
study. 

Geodetector is a statistical method used to detect spatial heterogeneity and its drivers. The core 
idea is based on the assumption that if an independent variable has a strong influence on a 
dependent variable, their spatial distributions should be similar (Wang et al., 2010). In this study, if 
an environmental factor influences glacier mass balance, the spatial distributions of the factor and 
the mass balance are expected to be similar. 

To probe the spatial stratification heterogeneity of variable Y, and the extent to which factor X 
explains the spatial heterogeneity of variable Y, we described the q value by Equation 2. 


L 
2 
YN, 


q = 1- h=1 (2) 
No? 


where N, and N are the unit number of strata h and the whole region, respectively; g? and o? are 
the variances of Y at strata h and the whole region, respectively; and the value range of q is 0-1. 
The q value implies that X explains 100g% of Y. The larger the q value, the more directly 
correlated X with Y. In other words, if the g value of the correlation between an environmental 
factor, and glacier mass balance is 1, then the factor's explanatory power is 100%. That is, the 
factor completely controls the glacier mass balance. 

A simple transformation of the description for the q value satisfies the non-central F 
distribution: 


peL EN, (3) 
L-1 1-q 
2 
ge |S yi Èv) (4) 
> h=- nth , 
o’ h=1 N kal 


where L is the maximum count of strata; ~ is the abeyance of some distribution; A is the non- 


central parameter; and Y, is the mean value of strata A. 


The principle of geodetector method is shown in Figure 1. The example diagram includes two 
coverages, dependent variable Y, and explanatory variable X. The explanatory variable X is a 
typical quantity expressed by a polygon. If X is grid data, it need be converted to type data by the 
reclassification method. In the classification, except when the slope is classified by eight 
directions, the classification scheme with the largest q value is followed. The dependent variable 
Y for every glacier is mass balance, and the average Y and variance (o?) are calculated. 


The two coverages are compared, enabling the calculation of the average Y value and variance 
(o) for each coverage. Substituting these mean values and variances into Equation 2 produces 
the g value. Subsequent substitution into Equations 3 and 4, and checking the non-central F 
distribution table results in the statistical significance test of the q value. 

Interaction between X; and X2 (gxinx2) is evaluated by comparing with values of qxı and qx, 
where gxinx2 represents the determinative force of the new factor resulting from the superposition 
of X; and X2 factors. All of the above processes are undertaken with the software of geodetector 
(www.geodetector.cn). 
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Fig. 1 Principle of geodetector. x, independent variable; y, dependent variable; h, the strata of x or y 
(classification or partition); Y, , the mean value of strata h; of, the variance of y at strata h; o°, the variance of y 


at the whole region. 


3 Results 


3.1 Glacier mass balance variability in the HMA 


The average glacier mass balance in the HMA was —0.18 (+0.04) m w.e/a during 2000-2016, but 
glacier mass change was heterogeneous (Fig. 2). For example, the glacier mass balance in the 
Hengduan and Nyainqêntanglha mountains was the most negative, which were —0.65 (+0.23) and 
—0.56 (+0.23) m w.e./a, respectively. The glaciers on the West Kunlun and East Kunlun mountains 
show a positive mass balance of 0.14 (+0.08) and 0.01 (+0.08) m w.e./a, respectively. The 
heterogeneity of the change of regional mass balance in the HMA is primarily correlated with 
climate. The glaciers in the southern part of the Tibetan Plateau have the most extensive mass loss 
due to the Indian monsoon. In contrast, the dominance of westerly wind results in the minimal 
glacier mass loss in northwestern HMA, and some glaciers are even positive. Glaciers in the 
central part of the Tibetan Plateau, dominated by transitional modes, had moderate glacier mass 
loss. The strengthening of westerly winds and the weakening of monsoon are the main drivers of 
the current glacial conditions in the HMA. Glacier morphology is also an essential factor affecting 
the heterogeneity of glacier mass balance. In particular, glaciers in the same region with similar 
climatic backgrounds show different change states. Two adjacent glaciers are identified as one 
having a positive balance and the other a negative balance. Therefore, in Section 3.2, we will 
discuss the influence of these environmental factors on glacier mass balance. 
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Fig. 2 Spatial distribution of the glacier mass balance change in the High Mountain Asia. Mts, Mountains. 
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3.2 Environmental factors affecting the variability of glacier mass balance 


Table 1 shows the influencing factors of the variability of glaciers mass balance in the HMA and 
other mountains, listing the q values from the largest to the smallest. The q value of DCR was the 
largest, inferring that the debris cover has an influential effect on the spatial differentiation of 
mass balance of the whole glaciers in the HMA. Moreover, climatic factors (T, Pc, and P) play a 
secondary role. Since not all glaciers are covered with supraglacial debris, climate should be the 
primary factor accounting for the spatial heterogeneity of glacier mass balance. However, the q 
value of temperature change (T.) is not higher, indicating that the spatial heterogeneity of 
temperature change has less influence on the spatial heterogeneity of glacier mass balance than 
other climate factors (T, Pe, and P). There is little difference in temperature change in the HMA, 
and spatial heterogeneity of temperature change has a relatively poor correlation with the spatial 


heterogeneity of glacier mass balance. 
Topography has less explanatory power for all the glaciers in the HMA than DCR and climate 
factors, but topography is the dominant effector of mass balance of most glaciers in mountain 


Region 


Table 1 Geodetector's results for different mountain regions 


Factors arranged from the largest to the smallest by the g value 


Dominant interaction (q) 


West Himalaya 
Central Himalaya 
East Himalaya 
Karakoram 


Pamir Plateau 
Hissar Alay 
Hindu Kush 
Western Tianshan 
Eastern Tianshan 


Western Kunlun 


Eastern Kunlun 


Tibetan Interior 
Mountains 


Tanggula 
Qilian 


Eastern Tibetan 
Plateau 


Gangdise Mountains 


Nyaingéntanglha 
Hengduan 


HMA 


S (0.20), NDG (0.08), Zmax (0.07), DG (0.06), P (0.05), Ps (0.03), T; (0.03), Zinea 
(0.03), PDG (0.03), T (0.03), Lmax (0.02), HI (0.02), A (0.02) 


NDG (0.18), S (0.18), DG (0.10), DCR (0.08), Zmea (0.08), HI (0.07), Zax (0.07), 
Zmin (0.07), Te (0.04), PDG (0.03), T (0.03), A (0.02), Lmax (0.02), P. (0.02) 


S (0.25), Zmea (0.17), Zmax (0.17), DCR (0.16), NDG (0.16), Zmin (0.15), P (0.13), 
DG (0.12), T (0.11), PDG (0.09) 


DCR (0.05), Zea (0.04), Zmax (0.04), T. (0.04), Zmin (0.04), P (0.03), T (0.03), S 
(0.03), Lmax (0.02), HI (0.02), A (0.02), NDG (0.01), P. (0.01), DG (0.01) 


S (0.13), Zmea (0.09), DCR (0.07), P. (0.06), NDG (0.06), A, (0.06), P (0.05), Zmax 
(0.05), DG (0.04), T (0.03), HI (0.03), T. (0.03), Zmin (0.03), PDG (0.02) 


P (0.16), HI (0.15), Te (0.15), S (0.14), T (0.12), Zmin (0.11), Pe (0.10), Lmax (0.09) 
S (0.41), DG (0.35), NDG (0.35), Zmax (0.24), PDG (0.17), Zea (0.13), P (0.11), 
P. (0.10), Zmin (0.09), Te (0.08), T (0.08), Lmax (0.07), HI (0.06) 


S (0.18), Zmea (0.17), Zmax (0.16), T (0.15), NDG (0.13), P (0.12), DG (0.11), Te 
(0.09), PDG (0.05), DCR (0.04), Zmin (0.04), Pe (0.04), A, (0.04), HI (0.03) 


Zmea (0.30), S (0.28), Zmax (0.27), NDG (0.25), DG (0.22), HI (0.15), Pc (0.11), T 
(0.08), Zmin (0.06) 


P (0.09), Emax (0.07), A (0.07), Zmea (0-05), DCR (0.04), Zmin (0.03), Te (0.03), Zmax 
(0.03), As (0.02), S (0.02) 


Zmed (0.28), Zmax (0.28), Zmin (0.26), T (0.23), Pe (0.23), P (0.15), DG (0.11), S 
(0.09), NDG (0.08), DCR (0.06) 


S (0.22), Zmea (0.21), Pe (0.21), P (0.20), Zmin (0.18), DCR (0.18), T (0.13), 
Lmax (0.09) 


Zmea (0.44), P (0.28), DCR (0.26), Zmin (0.25), Zmax (0.24), NDG (0.23), Te (0.22), 
T (0.18), P. (0.18) 


T. (0.17), P (0.15), T (0.14), Zmea (0.14), Zmax (0.08) 
Zmea (0.75), Zmax (0.72), PDG (0.71), DG (0.65) 


Zmed (0.52), Zmax (0.50), S (0.44), P. (0.43), Zmin (0.38) 


S (0.12), Lmax (0.10), P (0.08), T; (0.08), A (0.08), Zmea (0.07), T (0.05), Zmax (0.04), 
PDG (0.04), Zin (0.04), P. (0.04), A, (0.04), HI (0.03), DG (0.03), DCR (0.02) 
DG (0.49), Zmin (0.49), S (0.46), PDG (0.45), NDG (0.44), Zmax (0.38), Pe (0.29), 
T; (0.21), T (0.18), Zmea (0.18), HI (0.18), DCR (0.17), Lmax (0.16) 


DCR (0.14), T (0.10), P. (0.10), P (0.08), S (0.08), Zmea (0.06), Zmax (0.05), NDG 
(0.05), DG (0.03), HI (0.03), Te (0.02), Zmin (0.02), PDG (0.02), Lmax (0.01) 


SOP (0.25) 


NDGNDCR (0.33) 


SA Zinea (0-45) 


DCRNZmea (0.13) 


SNP. (0.22) 
TNT (0.41) 


SNNDG (0.55) 


ZmeaNT (0.28) 


SANHI (0.46) 


PA Lmax (0.18) 


Zmed! Pe (0.52) 


SMDCR (0.46) 


ZmeaM A (0.61) 
ZmeaNT (0.34) 
Zmeal IT (0.99) 
Zmin NT (0.79) 


SALmax (0.23) 


Zmin OAs (0.68) 


DCRNS (0.30) 


Note: A, area; S, slope; Zmin, minimum elevation; Zmax, Maximum elevation; Zmea, median elevation; As, aspect; Lmax, maximum length; DG, 
difference of glaciation; PDG, positive difference of glaciation; NDG, negative difference of glaciation; DCR, debris coverage rate; HI, 
hypsometric index; P, precipitation; T, temperature; Pe, precipitation change; T., temperature change; HMA, High Mountain Asia; We 
only keep statistically significant (at 99% level) the q values in the order from the largest to the smallest. Dominant interaction shows 
only the first pair of interactions arranged by the q values. 
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regions. Slope and median elevation have the highest g values among all topographic factors. 
Assuming that the climatic background of each mountain region is similar, slope and median 
elevation are the dominant factors of mass balance heterogeneity of glaciers in most mountain 
regions. However, there are some exceptions. For example, the slope has a weak correlation with 
the heterogeneity of glacier mass balance in the West Kunlun Mountains due to the generally 
gentle slope. Equally, the median elevation has a weak explanatory effect on the heterogeneity of 
glacier mass balance in the Hissar Alay Mountains due to the generally lower median elevation. 

Aspect, glacier area, maximum length and PDG have little influence on glacier mass balance, 
and the correlation is not significant (Tables 2 and 3). However, there are exceptions to this. For 
example, the glacier area, and the maximum length in the West Kunlun Mountains strongly 
correlate with the glacier mass balance. Glacier mass balance in the West Kunlun Mountains is 
positive, and the difference in glacier mass balance in this region is closely related to glacier size 
and precipitation. 

Our results show that the effects of any two environmental factors on glacier mass balance are 
mutually enhanced. Table 1 shows only the two environmental factors with the largest interaction 
enhancement effect (q value). For the whole HMA, the interaction between supraglacial debris, 
and the slope is the dominant. In combination with Tables 2 and 3, it is inferred that bare ice with 
no supraglacial debris coverage, and gentle slopes has a more negative mass balance. The pattern 
of the Tibetan Interior Mountains is the same as that in the whole HMA. The dominant interaction 
types of glacier mass balance in other mountain regions are different. However, in most cases, 
slope or median elevation is the most outstanding interaction compared with other factors. It is 
worth mentioning that the interaction between median elevations, and air temperature accounts 
for 99% of the glacier mass balance in the eastern Tibetan Plateau. 

The geographical analysis results (Table 1) of different mountain regions show that, except for 


Table 2 Correlation matrix analyzing glacier mass changes and environmental factors 


Region A Zmin Zmax Zmed S As Lmax DCR 
West Himalaya -0.11 -0.04 0.21 0.11” 0.45” 0.03 -0.11  —0.07 
Central Himalaya —0.09* —0.09* 0.22" 0.20" 0.43" 0.02 0.13"  —0.15" 
East Himalaya 0.00 0.03 0.34" 0.25" 0.48" -0.02 -0.03 -0.08 
Karakoram ~0.10"  —0.06* ~0.08"  —0.05* 0.14" 0.04 -0.13 013" 
Pamir Plateau ~0.01 0.08" 0.22” 0.25" 0.34" 0.11% -0.03 -0.06 
Hissar Alay —0.20* -0.13 -0.07 0.03 0.21" 0.06 -0.22" -0.17 
Hindu Kush 0.05 -0.19 0.43” 0.27" 0.63” 0.02 0.12 -0.09 
Western Tianshan —0.04 0.06 0.31" 0.40% 0.42™ —0.05 —0.05 0.05 
Eastern Tianshan 0.03 0.25"* 0.49" 0.54" 0.52” 0.08 0.07 0.12 
Western Kunlun —0.24** —0.10* -0.17" 0.14" 0.10" 0.16" 0.24" —0.10" 
Eastern Kunlun —0.13 0.39" 0.37" 0.53" 0.28" —0.06 —0.05 0.19"* 
beers -0.01 030" -0.16 -0.18 0.10 0.04 -0.02 -0.16 
Tanggula -0.02 0.19* 0.44" 0.64" 0.20" 0.03 -0.01 -0.31" 
Qilian -0.04 0.02 0.26" 0.28" 0.21* 0.05 0.02 -0.03 
E N Tibetan 0.18 -0.04 0.82" 0.85 017 0.24 0.30 -0.24 
Pac 0.01 0.52" 0.60" 0.70"* 0.59% 0.26 -0.16 -0.35" 
Nyaingéntanglha 0.24" 0.01 0.06 0.06 0.34" -0.01 -0.31" 0.02 
Hengduan Shan 0.26"  —0.58" 0.47" 0.02 0.63” 0.14 0.32" 0.32” 
HMA -0.00 0.11” 0.19" 0.20" 0.27” 0.01 -0.01 -0.19% 


Note: A, area; S, slope; Zmin, minimum elevation; Zmax, maximum elevation; Zmea, median elevation; As, aspect; Lmax, maximum length; 
DCR, debris coverage rate; HMA, High Mountain Asia; *, P<0.05 level; **, P<0.01 level. 
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Table 3 Correlation matrix analyzing glacier mass changes and environmental factors 


Region DG PDG NDG HI P T Pe Te 
West Himalaya 0.23"" 0.15” 0.25" -0.10% 0.03 0.07* 0.06 —0.09* 
Central Himalaya 0.26” 0.10* 0.38” -0.25% —0.03 0.06 0.04 0.00 
East Himalaya 0.29" 0.19" 0.31" —0.13* 021” 0.04 0.05 —0.07 
Karakoram —0.00 —0.03 0.03 —0.09** 0.05" —0.07™* —0.03 -0.11" 
Pamir Plateau 0.12" 0.05 0.17" -0.14" -0.13™" —0.12" -0.17" 0.11" 
Hissar Alay 0.04 —0.10 0.21* -0.37" -0.32™ —0.26"° —0.20° -0.29"" 
Hindu Kush 0.47" 0.35" 0.50” —0.20™ 0.24" 0.16" 0.25" 0.08 
Western Tianshan 0.22" 0.15” 0.27" -0.17" -0.21" -0.31" 0.11" —0.04 
Eastern Tianshan 0.40™ 0.17" 0.50 —0.33™ 0.02 —0.09 0.33 0.06 
Western Kunlun —0.04 —0.03 —0.04 0.01 -0.12" —0.08" —0.03 0.10" 
Eastern Kunlun 0.15" 0.07 0.20% —0.02 —0.10 —0.21™ —0.23™* 0.05 
Tibetan Interior Mountains 0.13 0.04 0.16 -0.09 -0.30" -0.17 —0.16 -0.11 
Tanggula 0.18" 0.14 0.18" —0.16 0.23" —0.20° 0.10 —0.00 
Qilian 0.17 0.08 0.23° —0.10 0.10 0.05 0.08 —0.20* 
Eastern Tibetan Plateau 0.65" 0.62" 0.56" 0.11 0.12 —0.24 0.22 0.06 
Gangdise Mountains 0.18 0.10 0.22 —0.06 0.14 —0.18 0.33* —0.32* 
Nyainqêntanglha —0.05 -0.14" 0.06 -0.15% 0.09* 0.06 0.14" 0.14" 
Hengduan 0.61” 0.55” 0.58” —0.11 0.04 0.20* 0.46” —0.28"" 
HMA 0.08"" 0.01 0.14"" -0.15% -0.23™" —0.29"* 0.10% -0.05™" 


Note: DG, difference of glaciation; PDG, positive difference of glaciation; NDG, negative difference of glaciation; HI, hypsometric index; 
P, precipitation; T, temperature; Pe, precipitation change; Te, temperature change; HMA, High Mountain Asia; *, P<0.05 level; *, 
P<0.01 level. 


the West Kunlun, Hissar Alay, and Qilian mountains, climate is not the dominant factor. This is 
presumably because the climate difference for each mountain region is slight, therefore, we 
assume that the climate background in each mountain region is similar. As a result of this 
assumption, the median elevation or slope becomes the dominant factor of glacial mass balance in 
regions with the same climatic background. The slope is the most crucial factor affecting glacier 
mass balance in the western Himalayas, eastern Himalayas, Hindu Kush, western Tianshan, 
Tibetan Interior Mountains, Nyainqéntanglha mountains, and Pamir Plateau. The median 
elevation is the most important factor affecting glacier mass balance in the eastern Tianshan, 
eastern Kunlun, Tanggula and Gangdise mountains, and eastern Tibetan Plateau. Precipitation 
plays a dominant role in glacier mass balance in the West Kunlun Mountains, where the 
correlation between the maximum length and precipitation is the strongest. Glaciers with longer 
lengths may receive more precipitation, increasing the influence of precipitation on glacier mass 
balance. The enhanced westerlies bring more precipitation to the West Kunlun Mountains, 
resulting in a positive glacier mass balance in the region (Yao et al., 2012; Lin et al., 2017). 
Precipitation plays a leading role in glacier mass balance in the Hissar Alay Mountains, but the 
correlation between temperature and temperature change is the strongest. Temperature change 
plays a dominant role in glacier mass balance in the Qilian Mountains, where the correlation 
between temperature and median altitude is the strongest. We also assessed the linear regression 
method, and the results are shown in Tables 2 and 3, which are also consistent with the results of 
Geodetector method. 


4 Discussion 


4.1 Limitations of this study 


The most significant limitation of this study is the accuracy of data, particularly the glacier mass 
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balance and climate reanalysis data. The individual glacier mass balance dataset is inherently 
highly uncertain, depending on the glacier area, the proportion of glacier surface in the survey, 
and the number of Digital Elevation Models (DEMs) used to extract reliable signal rate of 
elevation change (Brun et al., 2017). However, after the Monte Carlo simulations, Brun et al. 
(2019) suggested that the spatial heterogeneity of the data was relatively robust, and could be 
used for statistical analysis. Additionally, other data sources are also affected by the uncertainty of 
glacier delineation (Paul et al., 2013). Since glaciers larger than 2 km? were selected in this study, 
the glacier boundary uncertainty had little effect on the mass balance (Brun et al., 2019). The 
debris cover data, however, were uncertain, and the accuracy of the glacier classification was 91% 
(Kraaijenbrink et al., 2017), which was sufficient for this study. 

The spatial resolution (0.1°x0.1°) of climate reanalysis data is so coarse that it cannot reflect 
intra-mountain differences, limiting its accuracy. Although the value may be inaccurate, 
geodetector only need the classification of each glacier, not the exact value. We classify glaciers 
according to climate data, and divide them into eight types of climate conditions. Geodetector 
does not require accurate quantitative values, only qualitative classification results. The 
classification results can roughly reflect the differences in glacier climate. Our results (Table 1) 
show that climate in most mountainous regions has little influence on the difference in glacier 
mass balance, which is related to the coarse resolution of climate data. Thus, we can assume that 
climate condition in each mountain region is similar. Therefore, climate is the dominant factor 
affecting the difference of glacier mass balance in the HMA, which indicates that the resolution of 
climate data is sufficient to explain the difference in glacier mass balance in every mountain 
region of the HMA. In future studies, ERAS data could be scaled to the glacial level, and then 
divided into mountain regions to detect climate effects. However, using the new climate data does 
not change the results for other factors, because the q value of each factor in geodetector is 
calculated independently. The new calculations only affect the g values of climate factors which, 
in turn, affect the ranking of climate factors in relation to various environmental factors. 

Another limitation is that glacier surface features are not considered, such as a supraglacial 
lake, proglacial lake and ice cliffs. Ice cliffs and supraglacial lakes may exacerbate the glacier 
mass loss (Huang et al., 2018; Stefaniak et al., 2021), and these factors may complicate the effect 
of supraglacial debris. These effects are complex and need to be further explored in the future. 
Notwithstanding that if these factors will be considered in the future, then the impact factor q 
value of our present conclusions will remain the same. 

Finally, the current study fails to consider the temporal heterogeneity of glacier change. Studies 
have shown that the glacier mass in most regions shows an accelerated retreat trend, but the 
retreat degree is heterogeneous (Hugonnet et al., 2021). In terms of temporal heterogeneity, 
temperature and precipitation or their changes play essential roles, but topography and 
morphology cannot be ignored. 


4.2 Effect of supraglacial debris 


Due to the differences in reflectivity, particle size, color, and other physical properties, 
supraglacial debris has unique thermal characteristics that influence the physical processes of 
glacier ice ablation (Mattson et al., 1993; Nicholson and Benn, 2006). When the debris thickness 
is less than a critical threshold (<20—30 mm), the melting of underlying glacier ice is greater than 
that of bare ice; that is, the debris accelerates the melting of glacier. As the debris thickens (to 
greater than the critical threshold), the thermal resistance of the debris inhibits glacier ablation. 
Supraglacial debris is the strongest correlating factor explaining all the glacier mass balance of 
the HMA. However, many glaciers do not have supraglacial debris or are covered with a minimal 
amount of supraglacial debris. The reason is that the mean climate of a region is constrained by its 
geography, as well as its topography and morphology. There may have been insufficient sampling 
of debris-covered glaciers in various mountain regions, so supraglacial debris has not been 
reported to be the greatest influence on glacier mass balance in most mountain regions. Many 
glaciers are covered by supraglacial debris in the Karakoram Mountains; therefore, the 
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relationship between supraglacial debris coverage, and glacier mass balance has the highest 
reported q values in this region. And the supraglacial debris may be a contributing reason for the 
development of "Karakoram Anomaly". However, q values of all the influencing factors to the 
glaciers of the Karakoram Mountains are very low, which on the one hand is related to large 
numbers of debris-covered glaciers on the Karakoram Mountains, and on the other hand, may 
reveal the unique characteristics of the glaciers on the Karakoram Mountains (Table 1). 
Furthermore, the debris coverage also plays an essential role in explaining the glaciers mass 
balance heterogeneity of the Pamir Plateau, with the q value ranking behind the slope and 
median elevation factors (Table 1). 

Geodetector method does not easily differentiate whether the supraglacial debris has a positive 
or negative effect on the glacier mass balance. Therefore, a linear regression method was used to 
analyze this relationship. Except for the Nyaingéentanglha, Hengduan, and Tianshan mountains, 
the supraglacial debris of glaciers in other mountain regions are negatively correlated with the 
glacier mass balance (Table 2). However, only the Central Himalayas, Karakoram, East Kunlun, 
and Tanggula mountains are negatively correlated (P<0.01). The Hissar Alay, Western Kunlun, 
and Gangdise mountains are negatively correlated (P<0.05). The Hengduan Mountains is 
positively correlated (P<0.01). The thick supraglacial debris can inhibit the glacier mass balance, 
while the thin supraglacial debris can promote. The positive or negative correlation between 
supraglacial debris and mass balance may depend on the supraglacial debris thickness (Scherler et 
al., 2011; Benn et al., 2012; Fyffe et al., 2014). If the surface features of the debris are considered, 
such as ice cliffs and supraglacial lake, the influence of debris will be more complicated (Huang 
et al., 2018; Steiner et al., 2019; Stefaniak et al., 2021). Therefore, the effect of surface debris on 
mass balance is relevant and worth further exploration. 


4.3 Effect of climate 


Not all glaciers are covered with supraglacial debris, with many glaciers having minimal 
supraglacial debris. Therefore, climate is primarily dominant and significant factor for the 
heterogeneity of glacier mass balance throughout the HMA. Glacier mass balance is negatively 
correlated with temperature and precipitation (Table 2), indicating that glaciers in regions with 
higher temperatures or more precipitation may suffer a more severe mass loss under the 
background of climate warming. It seems counterintuitive that areas with much precipitation 
suffer more glacier mass loss. There are two reasons for this phenomenon. First, precipitation is 
negatively correlated with elevation, which means that areas with more precipitation tend to have 
lower elevations, and glaciers at lower elevations are more sensitive to rising temperatures 
(Oerlemans and Reichert, 2000). Secondly, there is more precipitation in the monsoon region and 
Indian monsoon-affected glaciers are extremely sensitive to temperature (Wang et al., 2019). The 
slight warming will lead to large-scale loss of glaciers, so Indian monsoon-affected glaciers are 
considered to be temperature-controlled glaciers. The effect of precipitation on glacier mass 
balance is complex. For the whole of HMA, the change of precipitation is positively correlated 
with glacier mass balance (Table 2). In particular, the increase of solid precipitation in recent 
years leads to a slightly positive balance of glaciers in the Karakoram Mountains (Kapnick et al., 
2014). In addition, the seasonality of precipitation varies greatly in the HMA, which affects the 
sensitivity of glacier mass balance s (Fujita, 2008). 


4.4 Effect of topography 


Topography may influence climate heterogeneity or climate sensitivity. Elevation influences 
temperature heterogeneity because it is warmer at lower elevations, and cooler at higher ones. Zmin, 
Zmax and Zmea are all related parameters that describe the elevation of glaciers. Aspect causes 
inhomogeneity in solar radiation, which may affect climatic heterogeneity (Anderson et al., 2010). 
Slope influences the flow and morphology of glaciers, thus, climate sensitivity. PDG, NDG, and 
HI are related to slope and elevation. Furthermore, the elevation range and slope determine the 
temperature and precipitation gradients on the glacier surface (Anderson and Mackintosh, 2012). 
Topography has less correlation for all the glaciers in the HMA than DCR and climate factors. 
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However, topography is the dominant factor for glacier mass balance in most mountain regions. 
Among all the topographic factors, median elevation, and slope are the most critical factors 
affecting glacier mass balance in each mountain region. The median elevation correlation with 
equilibrium line altitude and mass balance is zero at that elevation (Sakai et al., 2015). Hence, the 
median elevation is an important factor affecting glacier mass balance. The lower the median 
elevation, the lower the equilibrium line altitude, and the less the glacier mass loss. 

Slope-controlled glacier change has also been demonstrated by previous studies (Loibl et al., 
2014; Salerno et al., 2014). Data from current study shows that the gentler the slope, the more 
negative the glacier mass balance (Tables 2 and 3). Gentle slopes may lead to slower glacier flow, 
slower transport of glacial mass to the downstream, slower recharge in the ablation area, and more 
negative glacier mass balance. Moreover, the driving stress is conducive to separation and loss of 
ice (Salerno et al., 2014). 


4.5 Comparison of geodetector and linear regression analysis 


Univariate linear analysis was performed by calculating Pearson correlation coefficients, 
correlation of P values for mass balance and other environmental factors. The purpose of linear 
regression is to linearly and statistically relate the dependent variable Y to the independent 
variable X. Glacier mass balance is closely related to climate, environment, and terrain 
environment. The physical process of impact of climate and terrain on mass balance is very 
complex (Sakai et al., 2015). Their relationship is not necessarily linear, that is, not proportional, 
not linear quantitative relationship. Geodetector is suitable for linear and nonlinear relationship 
detection, so it is more suitable for exploring glacier mass balance factors than linear regression 
(Wang et al., 2010). And when the linear regression results are significant, the results of 
geodetector are inevitably significant. When the linear regression results are not significant, 
geodetector results may also be significant. Geodetector is more powerful, more confident, and 
strongly suggests causality better than general statistics. It is much more difficult for two 
variables to be uniformly distributed in two-dimensional space than for two variables to be 
uniformly distributed in one-dimensional curve. Therefore, we recommend that geodetector can 
be used in future studies to explore the influencing factors of heterogeneity in glacier mass 
balance. 


4.6 Sensitivity of glacier mass balance to environmental factors 


Sensitivity of glacier mass balance refers to the response degree of mass balance to temperature 
or precipitation change, which is generally the mass balance change when temperature increases 
by 1°C or precipitation increases by 10% (Oerlemans and Fortuin, 1992). Glacier mass balance is 
more sensitive to temperature at mid-to-high latitudes (Liu and Liu, 2016), and more sensitive to 
precipitation at lower elevations (Mölg and Hardy, 2007). This is because net shortwave radiation 
is the most variable energy flux at the glacier-atmosphere interface, and is controlled by surface 
albedo. Most of snow and ice mass loss (about 65%) occurs via sublimation (direct conversion to 
water vapor), followed by melting. The surface albedo depends on precipitation amount and 
frequency (Mélg and Hardy, 2007). Results has found that glacier mass balance has a nonlinear 
sensitivity to temperature and linear sensitivity to precipitation (Wang et al., 2019). Thus, 
temperature is the main driver of mass balance change in maritime glaciers, while in continental 
glaciers, the interactions among temperature, precipitation seasonality, and snow/rain separation 
are main drivers (Wang et al., 2019). Accordingly, topography may be the driver with obvious 
spatial heterogeneity, and nonlinear relationship with temperature. Thus, based on the temperature 
sensitivity classification of all glaciers in China by Wang et al. (2008), we used geodetector to 
explore the heterogeneity in temperature sensitivity affected by topography. The results show that 
the sequence of significant topographic factors is DCR (0.07)>Zmin (0.06)>Zmax (0.02)>HI 
(0.02)>Zmea (0.02)>slope (0.01)>aspect (0.01). Supraglacial debris can increase or decrease the 
rate of glacial melting, thus, affect temperature sensitivity (Mattson et al., 1993; Nicholson and 
Benn, 2006). The elevation range determines the temperature and precipitation gradients of 
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glacier surface (Anderson and Mackintosh, 2012). In this study, slope has less effect on 
temperature sensitivity than elevation. Slope sensitivity is not only controlled by temperature 
gradients but also by glacier flow and morphology (Loibl et al., 2014; Salerno et al., 2014). 
Glacier flow controls mass transport changes, which further lead to the response of changes in 
glacier mass balance. The glacier morphology affects the glacier flow. For example, when the 
glacier with a large accumulation area flows through a narrow valley bed, the flow velocity must 
increase (Salerno et al., 2014). Therefore, topography and morphology also control the spatial 
heterogeneity of mass balance by influencing the sensitivity of mass balance to climate change. In 
the future, more attention should be paid to the relationship among mass balance sensitivity to 
topographic and morphological parameters. 


5 Conclusions 


The purpose of this study is to examine the potential factors controlling the spatial distribution of 
glacier mass balance in the HMA during 2000-2016. The larger the q value of the influencing 
factor, the more significant the correlation of this factor to the glacier mass balance in the region. 
This research took advantage of various potential factors from available digital data, proposing a 
new geodetector computing framework to quantify, and compare glacier mass balance associated 
with various factors and the interaction between different factors. 

The results highlighted that supraglacial debris is an important factor affecting the spatial 
differentiation of glacier mass balance in the HMA, but not all glaciers are covered by 
supraglacial debris. Therefore, climate is the dominant factor for the spatial differentiation of 
glacier mass balance. When we assume that climatic background of each mountain is similar, the 
median elevation or slope is the dominant factor of the spatial differentiation of glacier mass 
balance. The effects of various factors on glacier mass balance are mutually enhanced, but there is 
no unified dominant and controlling interaction between factors. These outcomes suggest that 
geodetector provides a quantitative and objective analytical framework that can be used to 
investigate the potentially controlling factors of glacier mass balance more thoroughly. 
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